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ABSTRACT 

Under the assumption that the Milky Way's dark halo consists of primordial 
black hole MACHOs (PBHMACHOs), the mass density of the halo can be 
measured by the low frequency gravitational waves (10~ 3 Hz <; u gw ^ 1CT 1 Hz) 
from PBHMACHO binaries whose fraction is ~ 10~ 6 . We find that ten years 
observation by LISA will detect ~ 700 PBHMACHO binaries and enable us 
to determine the power index of the density profile within 10% (20%) and the 
core radius within 25% (50%) in about 90% (99%) confidence level, respectively. 
The axial ratios of the halo may also be determined within ~ 10%. LISA and 
OMEGA may give us an unique observational method to determine the density 
profile and the shape of the dark halo to open a new field of observational 
astronomy. 

Subject headings: black hole physics — dark matter — Galaxy: halo — Galaxy: 
structure — gravitation — gravitational lensing 
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1. INTRODUCTION 

It is important to determine the density profile of the Milky Way's dark halo 
observationally in order to gain insights into the galaxy formation and evolution. 
Unfortunately, little has been known about the halo density profile (HDP), since the dark 
halo emits little light. A HI rotation curve tells us about the radial distribution of dark 
matter (e.g. Fich & Tremaine 1991). However it has not been accurately measured for 
r > D ~ 8.5 kpc and recently there is also an argument that the Galactic rotation curve 
may deviate from that of the standard halo model (Honma & Sofue 1996, 1997, Oiling 
& Merrifield 1998). The dynamics of satellite galaxies and globular clusters can provide 
the mass inside r ~ 50 kpc with some biases (e.g. Lin, Jones, & Klemola 1995, Kochanek 
1996, Zaritsky et al. 1989, Einasto k Lynden-Bell 1982, Peebles 1995). The HDP may be 
recovered by using the tidal streams from Galactic satellites (e.g. Johnston et al. 1999). All 
methods so far are, however, indirect methods. In this paper we investigate a possibility of 
direct measurement of the density distribution of Galactic dark halo. 

Recently, Ioka, Tanaka, & Nakamura (1998b) (hereafter ITN) proposed a possibility 
to determine a map of a HDP by low frequency gravitational waves (10 _4 -10 _1 Hz) from 
PBHMACHO binaries, which can be detected by the planned Laser Interferometer Space 
Antenna (LISA) (Bender et al. 1998) and Orbiting Medium Explorer for Gravitational 
Astrophysics (OMEGA). ITN was motivated by the observations of gravitational 
microlensing toward the Large Magellanic Cloud (LMC). The analysis of the first 2.1 years 
of photometry of 8.5 x 10 6 stars in the LMC by the MACHO Collaboration (Alcock et 
al. 1997) suggests that the fraction 0.62io;2 of the halo consists of massive compact halo 
objects (MACHOs) of mass O.SIq^Mq assuming the standard spherical flat rotation halo 
model. At present, we do not know what MACHOs are. There are several candidates 
proposed to explain MACHOs, such as brown dwarfs, red dwarfs, white dwarfs and so on 
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(Chabrier 1999, Gould, Flynn, & Bahcall 1998, see also references in ITN). Any objects 
clustered somewhere between the LMC and the sun with the column density larger than 
25M pc~ 2 can also explain the data (Nakamura, Kan-ya, & Nishi 1996). They include the 
possibilities: LMC-LMC self-lensing, the thick disk, warps, tidal debris and so on (Sahu 
1994, Zhao 1998a,b, Evans et al. 1998, Gates et al. 1998, Bennett 1998, see also discussions 
on SMC events, Afonso et al. 1998, Albrow et al. 1998, Alcock et al. 1998, Honma 1999) 
However, none of them do not convincingly explain the microlensing events toward the 
LMC and SMC. 

Freese, Fields and Graff (1999) claimed that on theoretical grounds one is pushed to 
either exotic explanations or a non-MACHO halo. We here simply adopt the suggestion 
by the MACHO Collaboration (Alcock et al. 1997) and consider an example of exotic 
explanations: primordial black hole MACHO (Nakamura et al. 1997). This possibility is 
free from observational constraints at present (Fujita et al. 1998) and PBHMACHOs may 
be identified by LIGO, VIRGO, TAMA and GEO within next 5 years (Nakamura et al. 
1997, Ioka et al. 1998a), if they exist as dark matter. 

If primordial black holes (PBHs) were formed in the early universe at t ~ 1CT 5 s 
(Yokoyama 1997, Kawasaki & Yanagida 1999, Jedamzik 1997), a part of them evolved into 
binaries through the three body interactions (Nakamura et al. 1997, Ioka et al. 1998a). 
Some of these binaries emit gravitational waves (GWs) at low frequencies at present. ITN 
found that one year observation by LISA will be able to identify at least several hundreds 
of PBHMACHO binaries. Since LISA can measure distances and positions of PBHMACHO 
binaries (Bender et al. 1998, Cutler 1998), it may be possible to obtain HDP from the 
distribution map of PBHMACHO binaries. 

In this paper, we will quantitatively investigate how well HDP can be determined by 
the observation of the low frequency GWs from PBHMACHO binaries and show that LISA 
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and OMEGA will serve as excellent instruments for determination of the shape of our dark 
halo. 

2. PBHMACHO MODEL 

For simplicity, we assume that PBHs dominate the dark matter, i.e., Q = Qbh, where 
VLbh is the density parameter of PBHs at present, and that all PBHs have the same mass 
M B h- Throughout this paper, we will set M B h = O.5M and Vlh 2 = 0.1, where h is the 
present Hubble parameter H in unit of 100 km s _1 Mpc -1 . 

Assuming that PBHs are distributed randomly at their formation, we can obtain 
the probability distribution function (PDF) for orbital frequency v v and eccentricity e 
of the binary, f Vfi {y v ,e)dv v de (Nakamura et al. 1997, Ioka et al. 1998a). For e C 1, an 
approximate PDF is given by 

425 /to\ 3/37 /a\ 4 /58\ du p 

um-m*, p - — y ry-, a) 

where a = (GM B h /2vr 2 z/p) 1 / 3 is the semimajor axis, t = 10 10 years is the age of the 
universe, a = 2.0 x 1O 11 (M B ^/M ) 3//4 cm is the semimajor axis of a binary in a circular 
orbit which coalesces in f„, x = 1.2x 1O 16 (M B j / /M ) 1 / 3 (^/i 2 )- 4 / 3 cm is the mean separation 
of black holes at the time of matter-radiation equality and t = (3 7 (ax / a ) 4 £o (ITN). a and (3 
are constants of order unity. In this paper we adopt a = 0.5 and f3 — 0.7 (ITN). Note that 
a circular binary with orbital frequency v v emits GWs at GW frequency u gw = pu p with the 
second harmonic p = 2 (Peter & Mathew 1963, Hils 1991). 
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3. 



INDIVIDUALLY OBSERVABLE SOURCES 



To be observed as individual sources, the amplitudes of the GWs from the binaries 
have to exceed the threshold amplitude h th = 5/i J/ (Az/) 1 / 2 which is determined by the GW 
background h v and the frequency resolution Av = l/T (Schutz 1997, Thorne 1987). Here 
T is the observation time, and we set the signal-to-noise ratio (SNR) as 5.Q In our model, 
the GW background h u is determined by PBHMACHO binaries themselves (ITN, Hiscock 
1998). We use Figure 6 in ITN to estimate the GW background h v . The amplitude of the 
GW at v gw from a binary with eccentricity e, the harmonic p and the distance from the earth 
d is given by hi = 2^/GF i /c 3 7iv g ! w , where F { = L^(v gw , e)/4vrd 2 := L Q vfJ z p- lo l z g(p, e)/4vrd 2 
is the GW flux and L = (32c 5 /5G)(2ttGA^/c 3 ) 10 / 3 . The function g(p,e) is given by Peter 
& Mathew (1963), and M = M BH /2 1 ^ is the charp mass. Then, the requirement that the 
signal exceeds the threshold, hi > h th , determines the maximum distance to individually 
observable sources with p = 2 and e = as 



At frequencies above ~ v dis = 1(T 2 - 23 (T/lyear)" 6 / 11 (M BH /0.5M©)- 5 / 11 Hz, 
LISA can measure distances to PBHMACHO binaries since binaries change their GW 
frequencies v gw by more than Av through GW emission within the observation time 
T (Bender et al. 1998, Schutz 1986). Since binaries with high orbital frequencies 
(i/ p ^ 10~ 3 Hz) are almost in circular orbits at present (ITN), we can assume e = 
and p — 2. If a source with a circular orbit changes its GW frequency by £Az/ = £/T 

1 For simplicity we do not consider effects of the inclination of binaries and a reduction 
factor due to the antenna pattern of the detector in details. These effects can be absorbed 
in the observation time T, and the conclusion of this paper will hold if the effective T is 
increased correspondingly. 
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during the observation, the change in the binding energy of the binary is given by 
AE = (c 5 ^/3nu 2 gw GT)(7iu gw GM/c 3 ) 5 / 3 . From the relation L^(u gw , 0) = AE/T, we can 
obtain the charp mass as M = (57r^/96) 3 / 5 (c 3 /7rz/ 9U ,G)(7r^T)- 6 / 5 . Substituting the charp 
mass to the amplitude hi = (?>2/h) 1 l 2 (c/iiv gw d)('Ki> gw GM./ c^Y^ , we can obtain the distance 
to the source, d = (5/288) 1 / 2 (ct;/ir 2 b'g W h i T 2 ), as a function of hi, u gw , and £. 

The observable parameters, hi, v gm and £, contain observational errors as, hi(l ± ei), 
u gw (l ± e 2 ) and £(1 ± e 3 ), which can be estimated as e\ = 1/(SNR), e 2 = l/u gw T 
and £3 = l/£. Therefore the observational error in the distance can be estimated as 
d[l ± (ei + 3e 2 + e 3 )]. The angular resolution of LISA is estimated to be a few degree (Cutler 
1998)[[ 



4. DENSITY PROFILE RECONSTRUCTION 

In this section we show one of simulations of real observations. For simplicity, we 
assume that the distribution of the number density of PBHMACHOs obeys the law as 

= [! + («' (3) 
where r, n s , D a and A are the galactcentric radius, the number density of PBHMACHOs 
at the galactic center, the core radius and the power index, respectively. As the "real" 
parameters for a simulation we set n s = 2.60 x 10~ 2 pc~ 3 , D a = 5 kpc and A = 1. The 

2 The results of Cutler (1998) are only valid for large SNR and the angular resolution 
may be worse in a realistic detection with SNR = 5 (Balasubramanian, Sathyaprakash, & 
Dhurandhar 1996), although we here simply adopt the Cutler's results. Note also that the 
relative velocities of the sources to the solar system are not taken into account in Cutler 
(1998). 
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total number of PBHMACHOs within r < D ha i Q is given by N ta tai — I r <D halo n ( r )d 3 x 
where D^aio is the size of the halo . Since the fraction of PBHMACHO binaries with 
(1(T 3 Hz <^)v min < v gw < v max is given by F b ( 

Vmin i V m ax ) ■ 

with p = 2 from equation (|l|), the total number of PBHMACHO binaries with (10 -3 Hz 
<^)v m in < v gw < ^max is given by N b = N total F b (v min , u max ). Note that, as long as v gw ^ 10~ 3 
Hz, it is sufficient to consider the case e < 1 and the second harmonic p = 2. For example, 
N total = 4.03 x 10 12 , F b (u min} v max ) = 2.14 x 10" 8 and hence N b = 8.62 x 10 4 , for D halo = 500 
kpc, v min = 4 x 10 -3 Hz, and v max = 1 x 10" 1 Hz. 

The following algorithm explains a method of our simulations to see how well HDP can 
be determined by the observation of the low frequency GWs. 

1. We distribute N b PBHMACHO binaries randomly following the adopted HDP 
in r < D ha i and assigning frequencies according to the PDF in equation (|T|) for 

Vmin ^ V gw <C Vmax- 

2. We make an observation in this numerically generated galactic halo . Note that we 
cannot use all individual sources to reconstruct the density profile, since for low 
frequency a maximum distance to be observed as individual sources in equation (|2|) is 
short. If we want to determine the HDP within r < D obs , we have to use binaries with 
frequencies v gw > v obs where v obs is determined by D$ ax (v obs , 0) = D obs + D . Here 
Dq = 8.5 kpc is the distance from the galactic center to the earth. For simplicity, we 
use a uniform probability distribution to assign the observational error of e\ + 3e 2 + e 3 
to the distance and the error of 3° to the angular resolution. 

3. By fitting the distribution map of the HDP, we can compare the reconstructed HDP 
with the "real" HDP. Note that observationally the normalization n s in equation (|3]) 
has to be replaced by the density of PBHMACHO binaries at the galactic center n sb . 
n s is obtained from n s = n sb /F b (u obs , v max ). 
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An example of simulated observations is shown in Figure 1. This histogram shows 
the number iVj of the observed PBHMACHO binaries whose distances from the galactic 
center are in (i — l)5r < r < i5r (i = 1,2, ■ • •). We adopt 5r = 1 kpc to determine 
the structure within a few kpc from the galactic center. Here we set T = 10 years and 
D bs = 50 kpc, which corresponds to v ohs = 9.63 x 10 -3 Hz. In this realization, N map , the 
total number of PBHMACHO binaries which can be used to determine the HDP, is 719. In 
order to estimate the fitted parameters (n s &, D a and A), we apply the least squares methodP] 
minimizing x 2 = E* [A 7 * - ni(n s b, D a , X)] 2 /crf, where rii(n sb , D a , A) = A-Kn(r)dr. The 

variance <7j of A^ can be estimated by cr^ = ^/rij(n s {,, D a , A), since the distribution of iVj 
will follow the Poisson distribution with mean rii(n s b, D a , A) assuming that the statistical 
uncertainty dominates the instrumental uncertainty due to the observational errors. For 
this realization, the fitted parameters turn out to be n s /n r s eal = 0.779, D a = 6.20 kpc and 
A = 1.04, where n r s eal = 2.60 x 10~ 2 pc~ 3 is the "real" value. The reduced \ 2 is 0.913 with 
A7(=D b s /5r — 3) degrees of freedom. In Figure 2, the "real" parameter and the fitted 
parameter are marked with a filled square and a cross respectively in the D a -\ plane. The 
contours of constant A% 2 are also plotted with A% 2 = 1.00, 2.30, 4.00 and 6.17. In Figure 
3, the "real" HDP and the reconstructed HDP normalized by n r s eal are shown. It seems that 
in our method HDP is reconstructed quite well except for the central region. 

We performed 10 4 simulations of observations with (or without) the instrumental error 
to obtain the probability distributions of the core radius D a and the power index A. The 

3 Strictly speaking, we may have to maximize the probability for observing iVj 
PBHMACHO binaries in i-th bin from the Poisson distribution, P(n s b, D a , A) = 

^[ni(n s b, D a , \)] Ni e~ n '( n3b ' Da,x ' > / Ni\^ , instead of minimizing y 2 . However, since almost 
all Ni is larger than 10, it will be a reasonable assumption that the shape of the Poisson 
distributions governing the fluctuations is nearly Gaussian. 
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mean values (w) and the dispersions Aw = ((w 2 ) — (w) 2 ) 1 ^ 2 of these parameters w with 
(or without) the instrumental error are shown in Table 1. From Table 1, we find that 
the instrumental error does not affect the results so much.[] The probabilities that these 
parameters w are within \w — (w)\ < Aw and 2Aw turn out to be about 70% and 95% 
respectively from these realizations. Although the power index A is determined within 10% 
(20%) error in 89% (99.7%) confidence level (CL), respectively, by ten years observation, 
the dispersion of the core radius D a is somewhat large, 25% (50%) error in 63% (93%) CL. 

After we know the power index A accurately by this global observation, we can analyze 
the HDP for shorter distances r < D obs < D obs using the PBHMACHO bianries with lower 
frequencies u gw > obs , where obs (< u obs ) is determined by D<$ ax (P obs , 0) = D obs + D from 
equation (0). For example, for T = 10 years and D obs = 10 kpc, which corresponds to 
v bs = 5.15 x 10~ 3 Hz, the mean value and the dispersion of D a are found to be (D a ) = 4.81 
kpc and AD a = 0.710 kpc from 10 4 realizations with 5r = 0.5 kpc and A = 1. The 
dispersion AD a is reduced by a factor 0.5. Then, the core radius D a is determined within 
25% (50%) error in 91% (99.8%) CL. 



5. DISCUSSIONS 

In this paper we have quantitatively investigated how well the HDP consisting of 
PBHMACHOs can be determined by the observation of the low frequency GWs, assuming 
the spherical HDP in equation (Rt). We have found that ten years observation by LISA will 



4 This also justifies the assumption that the statistical error dominates the instrumental 
error. Note also that the dispersions Aw are consistent with the extent of the contour 
Ax 2 = 1-00 in Figure 2. This justifies the assumption that the distribution of the number 
Ni of the observed PBHMACHO binaries in i-th bin is nearly Gaussian. 
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be able to determine A, the power index of the HDP, within 10% (20%) error and D a , the 
core radius, within 25% (50%) error in about 90% (99%) CL, respectively. 

The halo of our galaxy may be non-spherical (e.g. Oiling and Merrifield 1997). For 
a non-spherical halo if we calculate quadrupole moments of positions of PBHMACHO 
binaries we can determine axial ratios ((c/a) and {b/a}) of the dark halo (for details see 
Dubinski & Carlberg 1991). Since about 700 PBHMACHO binaries can be used by ten 
years observation, errors in the axial ratios are estimated as less than 10% if the axial ratios 
are less than 0.8. 

We have assumed that MACHOs are PBHs. However they may be white dwarfs, or 
some other compact objects (e.g. Freese et al. 1999). For such cases also it may not be 
so strange to expect that some of them are binaries. If a fraction ~ 10~ 6 of them is in 
binary systems emitting GWs in the frequency range of 10~ 3 Hz ^ u gw <^ 10 _1 Hz, similar 
arguments to this paper will hold even for non-black hole MACHOs (see also Bond & Carr 
1984). 

We would like to thank professor H. Sato for continuous encouragement and useful 
discussions. We are also grateful to K. Taniguchi, T. Harada and H. Iguchi for useful 
discussions. This work was supported in part by Grant-in-Aid of Scientific Research of the 
Ministry of Education, Culture, Science and Sports, No. 9627 (KI), No.09640351 (TN), 
No.09NP0801 (TN). 
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Fig. 1. — The number iVj of the observed PBHMACHO binaries whose distances from the 
galactic center are within (i — l)5r < r < i5r (i — 1, 2, • • •) in one experimental realization 
is shown. We adopt 5r = 1 kpc. Here we set T = 10 years and D obs = 50 kpc, which 
corresponds to u obs = 9.63 x 10~ 3 Hz. The fitted curve (solid line) and the "real" curve 
(dashed line) are also shown. The fitted parameters are n s /n r s eal = 0.779, D a = 6.20 kpc and 
A = 1.04, where n r s eal = 2.60 x 10~ 2 pc~ 3 is the "real" value. The reduced x 2 is 0.913 with 
47(= D obs /5r — 3) degrees of freedom. 



-16- 



1.4 - 



1.2 - 



1 - 



n — i — | — i — i — i — | — i — i — i — | — 

T=10yr, D obs = 50kpc 

Ax 2 =1.00 

Ax 2 =2.30 

Ax 2 =4.00 

Ax 2 =6.17 



n — i — r 




X fit 
■ real 



0.8 



J I L 



J I L 



J I L 



J I L 



J I L 



6 8 
D a [kpc] 



10 



12 



Fig. 2. — The "real" parameter and the fitted parameter obtained from one experimental 
realization in Figure 1 are marked with a filled square and a cross respectively in the D a -\ 
plane. The contours of constant A% 2 are also plotted with A% 2 = 1-00, 2.30, 4.00 and 6.17. 
Note that for the Gaussian fluctuations the projections of the contours A% 2 = 1.00 and 4.00 
onto one axis contain 68.3% and 95.4% of data projected onto the axis respectively, and the 
contours A% 2 = 2.30 and 6.17 contain 68.3% and 95.4% of data respectively. 
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Fig. 3. — The "real" density profile with D a = 5 kpc and A = 1 (dashed line) and the density 
profile fitted from Figure 1 with n s /n r s eal = 0.779, D a = 6.20 kpc and A = 1.04 (solid line) 
are shown. These density profiles are normalized by n r s eal = 2.60 x 10~ 2 pc~ 3 . 
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Table 1. The mean values (w) and the dispersions Aw = ((w 2 ) — (w) 2 ) 1 ^ 2 of the core 
radius D a and the power index A obtained from 10 4 experimental realizations with (out of 
parentheses) and without (in parentheses) the instrumental error are shown for several 
observational time T. v^ s is the minimum frequency of the binaries to which LISA can 
measure distances. v ohs is the minimum frequency of the binaries which we can use to 
determine the density profile within r < D obs = 50 kpc. (N map ) is the mean number of the 
PBHMACHO binaries that can be used to determine the density profile. 



T 


Vdis 


Vobs 


(N ma p) 


(D a )±AD a 


(A) ± AA 


[year] 


[mHz] 


[mHz] 




[kpc] 




2 


4.08 


15.1 


217 (214) 


4.60 ± 2.44 (4.68 ± 2.43) 


1.00 ± 0.114 (1.01 ± 0.116) 


4 


2.80 


12.4 


367 (363) 


4.57 ± 1.90 (4.67 ± 1.90) 


0.992 ± 0.0876 (1.00 ± 0.0894) 


6 


2.24 


11.1 


496 (492) 


4.62 ± 1.64 (4.71 ± 1.64) 


0.992 ± 0.0763 (1.00 ± 0.0771) 


8 


1.92 


10.2 


614 (608) 


4.65 ± 1.47 (4.74 ± 1.46) 


0.991 ± 0.0684 (0.999 ± 0.0689) 


10 


1.70 


9.63 


723 (716) 


4.67 ± 1.35 (4.76 ± 1.34) 


0.990 ± 0.0626 (0.999 ± 0.0632) 



